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0.  INTRODUCTION 

Markov  processes  (with  discrete  state  space  over  continuous  time) 
are  good  models  for  many  stochastic  systems,  including  certain  queuing 
systems,  inventory  systems,  reliability  and  maintenance  models,  etc.  The 
analyses  of  these  Markov  process  models  are  often  restricted  to  steady- 
state  behavior,  that  is,  systems  in  equilibrium.  However,  there  are  many 
instances  in  which  the  transient  behavior  of  the  process  is  important: 

(i)  Transient  behavior  is  often  encountered  because  of  changes  in 

operating  conditions  of  the  system  due  to  exogenous  environmental 
conditions  or  internal  control  of  the  system,  such  as  the  opening 
or  closing  of  a  queuing  system. 

(ii)  In  some  systems  with  time  homogeneous  stochastic  behavior,  the 
convergence  to  steady  state  is  so  slow  that  the  equilibrium 
behavior  is  not  indicative  of  system  behavior. 
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(iii)  For  systems  in  steady  stale,  there  are  some  transients  of  interest, 
such  as  busy  periods  for  queuing  systems. 

(iv)  To  obtain  equilibrium  results  for  some  models,  transient  behavior 
is  useful.  For  example,  the  steady-state  behavior  of  regenerative 
processes  is  characterized  by  the  behavior  during  an  individual 
cycle  (transient  behavior);  see  Ross  [1970,  Theorem  5.8]. 

The  primary  transient  quantities  of  interest  are  the  time-dependent 
state  probabilities.  From  these,  time-dependent  moments  can  be  calculated 


etc.  (Other  quantities  of  interest  are  distributions  of  first  passage 


times,  expected  cumulative  occupancy  times  during  a  time  interval,  and 
expected  number  of  a  certain  class  of  events  occurring  during  a  time 
interval.)  In  general,  the  transient  state  probabilities  of  a  Markov 
process  can  be  computed  by  solving  a  system  of  linear  first  order  dif¬ 
ferential  equations.  Nice  analytical  solutions  rarely  exist.  Even  for 
the  M/M/1  queue,  the  solution  comes  out  as  a  rather  formidable  expression 
in  terms  of  modified  Bessel  functions;  see  Gross  and  Harris  [1974],  In 
general,  numerical  methods  must  be  used  to  solve  for  the  transient  state 
probabilities. 

This  paper  focuses  on  the  approach  called  "randomization,"  first 
introduced  as  a  computational  method  by  Grassmann  [1977].  It  is  a  general 
method  for  computing  transient  probabilities  of  Markov  processes  with 
finite  state  spaces,  for  example,  finite  queues  (finite  source  or  finite 
capacity).  The  method  can  be  used  on  many  infinite  state  Markov  processes 
by  approximating  the  process  with  a  finite  state  Markov  process;  for 
example,  an  M/M/c/00  queue  can  be  approximated  by  an  M/M/c/K  queue  with 
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sufficiently  large  K.  The  randomization  approach  has  the  advantage  of 
having  a  probabilistic  interpretation  that  can  be  exploited  in  efficient 
model  generation  and  computations  for  a  broad  class  of  Markov  models. 

The  paper  is  organized  as  follows.  Section  1  presents  transient 
theory  for  Markov  processes  and  reviews  some  of  the  solution  approaches. 
Section  2  presents  the  idea  of  subordinating  a  Markov  chain  to  a  Poisson 
process,  the  basis  of  the  randomization  approach.  Section  3  develops 
these  ideas  for  birth-death  processes.  Section  4  gives  details  encountered 
in  randomization  computational  algorithms.  Section  5  presents  a  general 
modeling  approach  called  5ERT  which  can  be  used  with  the  randomization 
approach;  three  examples  are  given.  Section  6  gives  formulae  and  meth¬ 
odology  for  computing  other  transient  quantities  in  addition  to  state 
probabilities.  Section  7  gives  a  randomization  implementation  for  sparse 
systems.  Section  8  contains  some  concluding  remarks. 

1.  TRANSIENT  STATE  PROBABILITIES  OF  MARKOV  PROCESSES 

Let  {X(t),  t  ^  0}  be  a  Markov  process  on  a  finite  state  space 
S  =  {0,1,..., N}  .  The  stochastic  evolution  of  this  process  is  character¬ 
ized  by  its  generator. 


where  the  rows  sum  to  zero,  i.e.. 
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and  the  element  q..  is 
1J 


y  q  .  .  -  q  .  =  o 

A-  1 


q..  =  lim  p..(At)/At  , 
1J  At-K)  1J 


where 


p.^At)  =  Pr{X(t  +  At)  =  j  |  X(t)  =  i}  , 

These  Markov  processes  include  realistic  classes  of  Markovian 
queuing  models  such  as  finite  waiting  room  (capacity)  queues,  finite 
source  (machine  repair)  queues,  cyclic  queues,  and  closed  queuing  net¬ 
works.  (Also,  infinite  Q  matrices  can  be  approximated  by  truncation  for 
some  large  value  of  N  for  many  models.)  For  example,  the  Q  matrix  for  a 


general  birth-death  system  with  finite  capacity  is  given  by 

0  0  0  ... 


...  o 


iix  (-VV 


0  0  . . . 


(  ^2  ^  • • * 


(‘Vi%i) 


Let  TTj(t)  =  Pr{X(t)  =  j}  and  Tr(t)  =  (TrQ(t),  tt^U),  ...,  tt^ ( t ) )  ; 
then  the  Kolmogorov  forward  equations  (see  Karlin  and  Taylor  [1975,  p.  152], 


or  Heyman  and  Sobel  [1982,  p.  295],  for  example)  are 


TT'(t)  =  Tr(t)Q  .  (2) 

For  a  finite  state  space  S  =  {0,1,..., N}  we  have  a  system  of  N  +  1 
first  order,  linear  differential  equations  to  solve. 

Given  ir(0)  ,  there  are  several  approaches  to  the  solution  of  (2). 
Classical  numerical  analysis  techniques  for  the  numerical  integration  of 
(2)  include  Runge-Kutta,  or  predictor-corrector  methods,  etc.,  as  dis¬ 
cussed  by  Arsham,  Balana,  and  Gross  [1981],  Grassmann  [1977],  or  Maron 
[1982].  These  methods  have  the  disadvantage  that  they  are  formal  numeri¬ 
cal  analysis  techniques  and  ignore  any  exploitable  probabilistic  structure 
of  the  problem  (for  example,  the  structure  of  Q). 

Another  approach  is  to  note  that 


7T ( t )  =  n(0)  e 


Qt 


is  the  solution  of  (2),  where  by  definition 


,Qt  _ 


-  1+  l 


*  Qntn 


n=l 


n! 


A  common  method  of  computing  e is  based  on  the  spectral  representation 
of  Q  (see  Jinlar  [1975,  pp.  367-370],  or  Karlin  and  Taylor  [1975,  pp. 
539-541],  for  example).  This  involves  computation  of  all  the  eigenvalues 
and  eigenvectors  of  Q,  a  very  tricky  feat  when  the  state  space  S  is  large. 
Liou  [1966]  proposed  an  iterative  computational  procedure  for  computing 


the  terms  of 


7r(t)  =  l  tt(0) 
n=0 


sV 

n! 


(3) 


recursively,  but  as  Grassmann  [1977]  points  out,  the  existence  of 
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negative  diagonal  elements  in  Q  leads  to  possibly  severe  round  off  errors. 
Moler  and  Van  Loan  [1978]  give  a  survey  of  methods  for  computing  the  expo¬ 
nential  of  a  matrix.  They  are  primarily  conceried  with  matrices  that  are 
small  enough  to  be  stored  in  core  memory  (the  randomization  method  can 
handle  much  larger  matrices  for  sparse  systems) .  They  also  point  out  how 
difficult  it  is  to  solve  this  computational  problem  in  general. 

The  randomization  method  is  especially  tailored  for  solution  of 
(2)  when  Q  is  the  generator  of  a  Markov  process  and  I  is  a  probability 
vector.  It  has  a  probabilistic  interpretation  that  can  be  used  to  bound 
the  truncation  error  arising  from  the  infinite  series  in  (3)  and  which 
allows  the  algorithm  to  work  only  with  positive  numbers  (thus  minimizing 
round  off  errors).  Furthermore,  the  probabilistic  meaning  is  an  aid  in 
setting  up  the  model  and  executing  the  algorithm  for  a  large  class  of 
systems.  The  randomization  approach  is  based  on  the  subordination  of  a 
Markov  chain  to  a  Poisson  process. 

2.  MARKOV  CHAINS  SUBORDINATED  TO  POISSON  PROCESSES 

Let  (Y^,  n  =  0,1,2,...}  be  a  Markov  chain  on  a  countable  state 

space  5  with  transition  matrix  P.  Let  (N(t),  t  0}  be  the  counting 

process  of  a  Poisson  process  with  rate  A,  i.e.,  N(t)  equals  the  number  of 

occurrences  in  [0,t]  and  Pr{N(t)  =  n}  =  e  ^t(At)n/n!  .  Assume  that  {Y  }  and 

n 

(N( • ) }  are  independent.  Define  the  process  {X(t)  =  t  >  0}  .  It 

is  well  known  that  (X(t),  t  >  0}  is  a  Markov  process  (continuous  time) 
on  S  with  generator  Q  =  A(P  -  I)  and  the  same  initial  distribution. 

This  construction  is  referred  to  as  a  Markov  chain  subordinated  to  a 
Poisson  process;  see  Feller  [1971,  p.  345],  or  £inlar  [1975,  p.  236]. 

-  6  - 
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Cohen  [1969]  refers  to  these  processes  as  derived  Markov  chains.  It  is 
also  occasionally  referred  to  as  a  Markon  chain  on  an  underlying  Poisson 
process.  Feller  [1971]  refers  to  the  events  of  the  Poisson  process  as  a 
"randomized  operational  time"  and  the  construction  as  "randomization  by 
Poisson  distributions;"  hence  the  name  randomization  is  applied  to  this 
construction  as  well  as  to  the  computational  technique  that  exploits  it. 

Conversely,  let  {X(t),  t  >_  0}  be  a  Markov  process  on  a  countable 

state  space  S  with  generator  Q.  Furthermore,  assume  that  X(*)  is 

uniformizable,  i.e.,  the  diagonal  elements  of  Q  are  uniformly  bounded,  and 

let  A  =  sup  q.'.  Then  there  exists  a  Markov  chain  (discrete  time) 
ie-S  1 

{Y^,  n  =  0,1,...}  on  5  with  transition  matrix  P  =  Q/A  +  I  and  a  Poisson 
process  (N(t),  t  0}  with  rate  A,  which  are  independent  of  each  other, 

such  that  the  process  t  _>  0}  and  (X(t) ,  t  >  0}  have  the  same 

finite  dimensional  distributions  and  are  thus  probabilistically  identical. 
(See  Heyman  and  Sobel  [1982,  pp.  310-311],  for  example.) 

The  above  construction  gives  a  very  useful  computational  formula 
for  the  transient  probabilities  of  a  uniformizable  Markov  process.  Con¬ 
ditioning  over  the  number  of  occurrences  of  the  Poisson  process  in  [0,t] 
and  using  the  law  of  total  probability  gives 

tt  (t)  =  P{X(t)  =  s} 
s 

■  P{YN(t)  ■  S> 

<x> 

=  l  p(YN(t)  =  s  )  N (t)  =  n}  P{N(t)  =  n} 


■  I  p,v„  '  s) 

11=0 


e-At(At)n 
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Defining  4>s(n)  =  *’^Yn  =  and  cf> (n)  =  (<{>^(n).  ({^(n),  ...)  ,  this  equa¬ 
tion  becomes 

.<«> .  'f  *<„>  dhpl .  (4) 

n=0  ~ 

Recalling  that  n(0)  =  <J>(0)  and  P  is  the  transition  matrix  of 
(Y^,  n  =  0,1,...}  ,  equation  (4)  can  be  written  as 

u(t)  =  l  it (0) pn  - - .  (5) 

n=0  ~  n* 

This  formula  appears  in  Cinlar  [1975,  p.  259],  and  Cohen  [1969,  pp.  46- 
47].  It  is  apparently  originally  due  to  Jensen  [1953]. 

This  construct  relating  a  Markov  process  to  a  randomized  Markov 
chain,  originally  due  to  Jensen  [1953],  is  useful  in  a  least  two  other 
general  areas  of  applied  probability.  It  is  useful  for  developing  sto¬ 
chastic  inequalities  between  processes  and  monotonicity  properties  of 
processes  (see  Keilson  [1979],  Keilson  and  Kester  [1977],  Kirsten  [1977], 
and  Sonderman  [1980],  for  example).  This  construction  is  also  useful  for 
establishing  equivalences  between  discrete  and  continuous  time  Markov 
decision  processes  (see  Howard  [1960],  Lippman  [1975],  Serfozo  [1979], 
and  Veinott  [1969]). 

3.  A  DIRECT  PROBABILISTIC  DEVELOPMENT  OF  THE  RANDOMIZATION 
CONSTRUCT  FOR  A  BIRTH-DEATH  PROCESS 

As  an  example  to  illustrate  this  construction,  let  us  consider 
the  birth  and  death  process  whose  generator  Q  is  given  in  (1).  We  have 
a  birth-death  process,  which  means  that  the  system  state  (the  population 
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size)  can  change  by  +1  or  -1;  that  is,  if  a  birth  occurs,  the  size  in¬ 
creases  by  one,  while  if  a  death  occurs,  the  sis:e  decreases  by  one.  In  a 
small  instant  of  time,  if  the  system  is  in  state  s,  the  probability  of  an 
increase  to  s  +  1  is  AgAt  +  o(At)  ,  and  the  probability  of  a  decrease  to 
s  -  1  is  ygAt  +  o(At)  .  By  the  definition  of  conditional  probability, 
letting  X(t)  represent  the  random  variable  "size  of  the  population  at 
time  t,"  we  have 

Pr{X(t  +  At)  =  s  +  1  |  X(t)  =  s  and  a  transition  occurs  in  (t  +  At]} 

[  A  At  +  o(At)  ] it  (t) 

_ s _ _ _ s _ 

[A  At  +  y  At  +  o (At )  ] tt  (t) 
s  s  s 

A  At  +  o(At) 

_ s _ 

(A  +  y  )At  +  o(At) 
s  s 


Taking  the  limit  as  At  goes  to  zero  tells  us  that  given  that  a  transition 

occurs  when  the  system  is  in  state  s,  the  probability  that  it  will  be  a 

birth  (increase  to  s  +  1)  is  A  /(A  +  y  )  ,  where  y_  =  0  and  A„  =  0  . 

s  S  S  (J  N 

Thus  the  probability  that,  given  a  transition  occurs  when  the  system  is 

in  state  s,  it  is  a  death  (decrease  to  s  -  1)  is  y  (A  +  y  )  . 

s  s  s 

The  holding  time  in  state  s  for  this  Markov  process  is  the  minimum 
of  two  exponential  random  variables,  the  interarrival  time  with  parameter 


A  and  the  service  time  with  parameter  y  ;  hence  the  holding  time  is  also 
s  s 

an  exponential  random  variable  with  parameter  Ag  +  yg  . 

Thus  we  can  view  the  system  as  a  Markov  process  with  exponential 

holding  times  [mean  of  (A  +  y  )  ^  for  holding  in  state  s]  and  transition 

s  s 

probabilities  of  A  /(A  +  y  )  and  y  /(A  +  y  )  for  going  "up"  one  or 

s  s  s  s  s  s 

"down"  one,  respectively,  from  state  s.  Furthermore,  it  can  be  shown 
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that  the  occurrence  of  "birth"  or  "death"  is  independent  of  the  holding 
time. 


Suppose  we  wish  to  generate  this  process,  for  example,  to  simulate 
it  by  Monte  Carlo  means.  We  could  first  generate  the  transition  time  oc¬ 
currence  and  then  generate,  using  simple  Bernoulli  probabilities,  the 
change  of  state — that  is,  whether  the  process  increases  or  decreases  its 
state  by  one.  The  generation  of  the  transition  times  is  somewhat  compli¬ 
cated  in  that  we  would  have  to  sample  from  an  exponential  distribution 

with  a  state-dependent  parameter  X  +  u 

s  s 

To  avoid  this,  we  can  reproduce  this  process  in  the  following  way. 

Consider  the  minimum  mean  holding  time  (time  until  the  next  transition 

occurs)  of  the  process.  This  will  correspond  to  the  maximum  value  of 

A  +  p  ,  or  equivalently,  the  minimum  diagonal  element  of  Q,  the  gener- 
«s  s 

ator  of  the  process.  Call  this  value  A-  Let  us  denote  the  diagonal  ele¬ 
ments  of  the  Q  matrix  by  -qg,  that  is. 


qs  = 


X  +  p 
s  "s 


lMN 


(s  =  0) 

(s  =  1,2, . . . ,N  -  1)  . 
(s  =  N) 


To  reproduce  our  desired  transition  occurrence  process,  we  could  generate 
a  true  Poisson  process  with  rate  A  (exponential  holding  times  with  con¬ 
stant  mean  1/A)  and  then  "thin"  the  process  to  get  the  desired  state- 
dependent  transition  rate.  By  thinning,  we  mean  that  whenever  an  occur¬ 
rence  is  generated  by  the  Poisson  (A)  process,  if  we  are  in  state  s  we 
draw  from  a  Bernoulli  probability  distribution  with  success  probability 
qg/A,  where  a  success  indicates  counting  the  occurrence  of  a  transition 
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and  a  failure  (probability  1  -  qg/A)  indicates  ignoring  the  occurrence. 

This  thinning  procedure  on  the  Poisson  (A)  process  generates  our  desired 

underlying  state-dependent  transition  process. 

Thus,  to  reproduce  our  process  we  can 

i.  generate  a  Poisson  process  with  rate  A; 

ii.  thin  the  Poisson  (A)  process  by  a  Bernoulli  "switch"  with 

acceptance  probability  q  /A  and  rejection  probability 

s 

1  -  qs/A; 

iii.  generate  the  state  change  (up  or  down  one  unit)  by  another 

Bernoulli  switch  with  an  up  probability  of  X  / (A  +  y  ) 

s  s  s 

and  a  down  probability  of  y  / (X  +  y  )  . 

s  s  s 

Denote  the  Poisson  process  as  (N(t),  t  >  0}  ,  i.e.,  N(t)  equals  the 
number  of  occurrences  in  [0,t],  and  let  equal  the  state  of  the  system 
after  the  nth  occurrence  of  the  Poisson  process,  i.e.,  Y^  equals  X(Tn>, 
where  T^  =  min(t  :  N(t)  n)  .  Viewing  the  process  in  this  manner  will 
allow  us  to  derive  the  randomization  computing  algorithm  that  will  yield 
TT(t). 

Consider  an  element  of  the  transient  state  probability  vector 
rr(t),  say  ng(t)  =  P{X(t)  =  s}  .  Letting  pig(t)  =  Pr{X(t)  ■  s  |  X(0)  -  i} 
gives 

"  (t)  -  l  7r  (0)p  (t)  .  (6) 

s  ieS  s 

Now  considering  our  process  as  described  above,  we  have,  using  the  law  of 


total  probability. 
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F', 

i  • 


Pis(t)  =  2  Pr{Yn  =  8  I  Yo  “  i}  =  n} 

n=0 


00  ,  v  -At,,  .n 

_  £  p(n)  e  (At) 


n=0 


is 


n! 


(7) 


~(n) 


where  p.  represents  the  probability  of  going  from  state  i  to  state  s 

IS 

in  n  occurrences  of  the  Poisson  process.  But 

p^^  =  Pr{of  going  from  i  to  s  in  one  occurrence} 

=  Pr{of  accepting  the  occurrence  as  a  transition} 
x  Pr{of  transiting  from  i  to  s  I  a  transition  takes  place} 


1  - 

Vo 


xi  +  pi 


K  +  Va 


Ai  +  Pi 
A 


~r  (s  -  i  +  i) 


«  IT  (s  =  i  -  1) 


(s  =  i) 
(elsewhere) 


cU> 


Denote  the  matrix  with  elements  p^g'  by  P;  note  that  P  =  Q/A  +1.  We 
can  obtain  as  elements  of  a  matrix  formed  by  multiplying  P  by  itself 

n  times,  i.e. , 

ip[f}  =  [P]n  . 

Substituting  this  into  equations  (6)  and  (7)  yields 

■At . 


w(t)  =  l  tt(0)P 
0=0"” 


~n  e 


(At) 


n 


n! 


which  is  equation  (5) .  Thus  we  have  illustrated  the  general  theory  with 
this  special  case  of  a  birth-death  process. 
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4.  ALGORITHMIC  IMPLEMENT AT ION  OF  THE 
RANDOMIZATION  FORMULA 


The  randomization  algorithm  is  based  upon  equation  (4): 

n(t)  =  l  <t»(n)  .  (4) 

n=0 

There  are  several  considerations  that  arise  in  implementing  this  tech¬ 
nique:  (i)  truncation  of  the  infinite  series,  (ii)  efficient  calculation 
of  <j>n,  (iii)  multiple  time  points,  (iv)  changes  in  Q,  (v)  an  approach  to 
systems  with  small  state  spaces,  (vi)  an  approach  to  systems  with  infi¬ 
nite  state  spaces,  and  (vii)  efficient  computation  of  moments  and  prob¬ 
abilities  of  subsets. 


4.1  Truncation  of  the  Randomization  Formula 

Equation  (4)  involves  an  infinite  series.  Thus,  we  must  obviously 
truncate  the  series  at  some  point,  say  T,  so  that  the  computational  formu¬ 
la  becomes 


TT(t)  =  I  <J>(n) 
n=0 


-At,.  >n 
e  (At) 


nf 


(8) 


We  can  set  T  to  bound  the  error  due  to  truncation  (see  Arsham,et 
al.  [1981],  for  example)  by  satisfying 

T  IMn 

n  I 

n=0 


-  ^  I  e  , 

_ A  41  • 


where  e  is  the  desired  error  control.  This  would  guarantee  it  (t)  and  all 

s 


.-3 


other  probabilities  to  be  accurate  to  within  e  (e  *  10  ,  for  example, 

guarantees  at  least  two  decimal  place  accuracy) . 
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4.2  Efficient  Calculation  of  <j)  (n) 

To  use  equation  (8),  4> (n)  must  be  computed.  Recall  that  <J)g (n)  = 

P{Y  =  s},  where  {Y  ,  n  *  0,1,...}  is  the  discrete  time  Markov  chain, 
n  n 

Thus  the  randomization  computational  technique  really  reduces  to  compu¬ 
tation  of  the  transient  probabilities  <J>  of  the  Markov  chain. 

Computation  of  the  <j>(n)'s  is  naturally  a  recursive  procedure,  since 

<J>(n  +  1)  =  <fr(n)P  ,  (9) 

a  fact  that  follows  from  elementary  Markov  chain  theory.  In  general, 
the  calculation  in  equation  (9)  is  a  straightforward  matrix  multiplica¬ 
tion.  However,  most  specific  systems  have  sparse  P  matrices  and  thus 
efficient  computation  of  equation  (9)  should  exploit  the  sparseness  of 
P.  Grassmann  [1977,1982]  uses  sparse  matrix  techniques,  although  he 
does  not  explain  them  in  detail.  Melamed  and  Yadin  [1981]  give  a  recur¬ 
sion  for  efficiently  computing  4> (n)  for  Jacksonian  queuing  networks  which 
exploits  the  sparseness  of  Q.  We  give  a  general  modeling  approach 
called  -SERT  in  Section  5  of  this  paper — it  is  a  method  applicable  to  many 
sparse  systems.  In  Section  7  we  give  a  general  approach  that  completely 
exploits  the  sparseness. 

There  is  one  other  feature  of  the  efficient  calculation  of  equa¬ 
tion  (9)  that  should  be  mentioned.  In  many  systems,  it  is  most  natural 
to  use  a  multidimensional  description  of  the  state  space.  However, 
equation  (9)  can  be  calculated  in  less  time  if  the  state  space  is  "lin¬ 
earized,"  i.e.,  the  states  are  ordered  into  a  one-dimensional  state  space. 
This  efficiency  is  gained  because  less  arithmetic  is  involved  in  locating 
and  storing  values  in  a  vector  than  in  an  array.  Also,  since  these 
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multidimensional  state  spaces  often  have  irregular  shapes,  some  storage 
space  is  wasted  by  the  dummy  states  required  in  filling  the  state  space 
out  as  a  rectangular  array. 

4.3  Transients  at  Multiple  Time  Points 

In  a  typical  application  we  will  want  to  compute  transient  quan¬ 
tities  at  a  sequence  of  time  points  t^  <  <  ...  <  <  t^  <  t^ ^  <  t^ 

instead  of  just  at  a  single  time  point.  Typically,  these  time  points  will 

be  on  a  lattice,  i.e.,  t  -  t  *  h  for  all  m.  Without  loss  of  gener- 

m+l  m 

ality,  we  can  label  our  time  scale  so  that  h  *  1.  Consequently,  we  would 
like  to  compute  7r(l),  tt(2),  ...,  ff(£)  and  related  measures  over  the 
time  periods. 

Using  the  randomization  algorithm  Z  times,  once  for  each  ir(m),  is 
very  inefficient  because  the  <j)'s  will  be  recomputed  each  time.  Thus, 
the  efficient  approach  is  to  compute  the  <J>'s  recursively  and  as  each  is 
computed,  add  a  term  to  each  of  the  Z  series 

-At 

1  e  m(At  )n 

I  <J>(n)  - -  -m--  -  ,  m  *  1,2, ...  ,Z  , 

n=0 

until  the  desired  truncation  points  T(e,t  ),  m  =  1,2, ... ,Z,  are  reached. 

m 

4.4  Changes  in  the  Generator  Q 

A  phenomenon  that  may  occur  in  the  study  of  transient  behavior  is 
time-nonhomogeneity  of  the  stochastic  process.  For  a  Markov  process,  one 
manifestation  of  this  phenomenon  is  the  change  of  the  generator  at  some 
discrete  point  in  time:  on  [0,£),  the  process  may  have  generator  Q  and 
then  on  [£,,£']  have  generator  Q'.  In  order  to  compute  the  transient 
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probabilities  in  this  case,  we  compute  it ( 1 ) . tt(£)  as  before  and 

then  use  rr(£)  as  the  initial  probability  vector  for  a  second  problem 
involving  a  process  with  generator  Q'. 

Thus,  letting  P'  *  Q'/A'  +  I  be  the  transition  matrix  of  a 

Markov  chain  {Y1,  n  =  0,1,2,...}  ,  which  has  initial  distribution 
n 

<J>'(0)  -  tt (H)  and  marginals  d>'(n)  =  Pr{Y'  =  s}  ,  we  have 
~  ~  s  n 

$'(n  +  1)  =  ♦'(n)P' 

and 

T'  -A*k..,..n 

ir(4  +  k)  =  l  <j)'  (n)<- - \  =*-  ,  k  =  1,2,...  . 

n=0  ~  n‘ 

When  applying  the  randomization  algorithm  to  Markov  processes 
whose  generator  changes  at  discrete  times  during  the  time  interval  of 
interest,  [0,?.],  there  is  a  slight  difficulty  in  dealing  with  the  trun¬ 
cation  error  e.  It  is  desired  that  all  probabilities  computed  have  a 
maximum  error  equal  to  the  inputted  value  of  e;  this  is  achieved  by 
working  over  each  subinterval  with  epsilons  which  sum  to  an  inputted  E 
and  are  proportional  to  the  appropriate  time  subinterval.  For  example, 
suppose  we  are  interested  in  transient  behavior  over  [0,15]  and  the 
generator  changes  at  times  6  and  10.  Then  the  problem  will  be  solved 
in  three  parts:  first  on  [0,6],  epsilon  equal  to  (6/15)e  is  used; 
then  on  [6,10],  epsilon  equal  to  (4/15)e  is  used;  and  finally  on  [10,15], 
epsilon  equal  to  (5/15) £  is  used.  This  approach  will  guarantee  that  all 
probabilities  related  to  the  interval  [0,15]  will  have  a  computational 
error  of  e  or  less. 
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4.5  Systems  with  Small  State  Spaces 


For  systems  with  small  state  spaces  which  do  not  have  a  high 

degree  of  sparseness,  there  is  another  way  to  compute  transient  prob 

abilities  on  a  lattice  of  time  points,  t  =  m,  m  =  0,1,..., i  .  In 

m 

this  case  we  modified  the  above  algorithm  to  compute  the  matrix  of 
probabilities  P(l)  =  {p^g(l)}  »  where 

p  (1)  =  P{X(1)  =  s  |  X(0)  =  i} 

1»S 


and  using  previous  notation  we  have 


°°  -A  .n 

P(l)  =  l  Pn  6  A 


n=0 


n ! 


giving  us  a  computational  formula  after  the  appropriate  truncation  of 
the  infinite  sum. 

Assuming  the  process  is  time-homogeneous  over  [0,Jt)»  i.e.,  the 
generator  Q  does  not  change,  we  use  the  initial  probability  vector  ir(0) 
and  then  recursive  computation  to  obtain 

r(l)  =  tt(0)P(1) 

7T(2)  =  TT(1)P(1) 


tt(JI)  =  tt( £  -  l)P(l)  . 

The  program  SASTRANX  described  in  Balana,  et  al.  [1982]  computes  tran¬ 
sient  probabilities  for  a  machine  repair  system  that  experiences  changes 
in  rates  at  a  finite  number  of  discrete  time  points  in  this  way. 
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! 

! 


J 

i 

to 


4.6  Systems  with  Infinite  State  Spaces 

Some  systems,  such  as  infinite  capacity  queuing  systems,  will 
have  infinite  state  spaces.  Melamed  and  Yadin  [1981]  consider  this 
type  of  system.  In  this  case  it  is  necessary  to  truncate  the  state 
space.  This  will  introduce  an  approximation  error  into  the  calcula¬ 
tions.  This  error  can  be  bounded  by  collapsing  all  the  truncated  states 
into  one  absorbing  state  and  then  applying  the  randomization  algorithm. 
The  probability  that  the  system  is  in  the  absorbing  state  will  be  a 
bound  on  the  error  introduced  by  truncating  the  state  space.  For  an 
example,  see  Section  5.2 

4.7  Efficient  Computation  of  Moments  and 
Subset  Probabilities 

In  the  transient  analysis  of  a  Markov  process,  if  the  only  per¬ 
formance  measures  of  interest  are  moments  or  probabilities  of  certain 
subsets  of  the  state  space,  then  it  is  unnecessary  to  compute  the  entire 
state  probability  vector  tt( t )  . 

Suppose,  for  example,  we  are  only  interested  in  P{X(t)  =  A)  for 
some  subset  A  of  the  state  space.  It  will  be  more  efficient  to  forego 
computation  of  ir(t),s  and  to  compute  probabilities  at  occurrence  times 
of  the  underlying  Poisson  process 

An  =  l  <t>s(n)  «  P(Xn  e  A}  =  P{X(Tn)  e  A) 
seA 

and  then  use 

T(e,t) 

A(t)  =  P(X(t)  e  A}  =  l  A  P{N(t)  =  n}  . 

n=0  n 
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Similarly,  transient  measures  of  the  form 


E[f  (X(t  )  1  m  =  1 , 2 ,....,  H 

m 

can  be  more  efficiently  computed  by  foregoing  ir(t  ),  m  = 

~  m 

instead  using 


and 


E [f (Y  )  1  =  l  f(s)4>  (n) 
n  seS  s 

T(e,t) 

E [f (X(t) ) ]  =  l  E[f(Y  ) ]P(N(t)  =  n)  . 
n=0 


1, 


. ,  l  and 


An  example  of  this  type  of  performance  measure  is  the  expected  number 
of  machines  operating  at  time  t  in  a  machine  repair  model. 


5.  SERT  MODELING  PROCESS  AND  RANDOMIZATION 
We  now  present  our  approach  to  modeling  continuous  parameter 
Markov  processes,  which  utilizes  the  randomization  procedure  in  an  effi¬ 
cient  manner  and  enables  us  to  determine  the  transient  quantities  of 
interest  discussed  in  the  previous  sections.  We  call  the  approach 
SERT,  since  it  involves  the  four  concepts, 

5  =  state  space 
E  =  set  of  types  of  events 
R  =  set  of  transition  rate  vectors 
(one  for  each  element  in  E;  each 
vector  has  |S|  components) 

T  =  set  of  target  state  vectors 

(one  for  each  element  in  E;  each 
vector  has  |S|  components) 
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We  now  describe  an  algorithm  based  on  the  above  for  a  general 

class  of  continuous- time  Markov  processes  on  discrete  state  spaces. 

Consider  a  Markov  process  on  a  discrete  state  space  S.  It  has  E 

types  of  events  that  may  occur:  e. ,  e„,  ....  e  .  Let  E  =  {e, ,e_, . . . ,e  } 

12  E  I  2  E 

For  each  e.  e  E  there  exist  two  vectors,  the  vector  of  transition  rates 
J 

r J  ,  and  the  vector  of  target  states  t  ,  j  =  1,2,...,E.  Thus,  when 
the  process  is  in  state  s  e  S,  r^,  the  sth  component  of  r J ,  is  the  rate 
at  which  event  e^  will  occur  and  consequently  put  the  system  into  state  t 
The  characterization  of  the  process  in  terms  of  states,  events, 
rates,  and  targets  is  a  very  fruitful  way  of  working  with  a  broad  class 
of  Markov  processes,  both  from  the  point  of  view  of  defining  or  describ¬ 
ing  the  process  and  for  use  in  algorithmic  solution  methodologies. 

The  holding  time  in  state  s  is  exponentially  distributed  with  rate 
E  i 

r  =  l-  ,  r.  Let  A  =  max  r  be  the  rate  of  an  underlying  Poisson  pro- 
S  J  scS  S 

cess  (this  corresponds  to  the  maximum  diagonal  element  of  the  Q  matrix 
but  with  this  approach  it  is  not  necessary  to  generate  the  Q  matrix  per 
se) .  Define  a  Markov  chain  on  this  Poisson  process  whose  transitions 
are  the  result  of  E  +  1  distinct  events;  e^  is  the  "null  event."  The 
transition  probabilities  corresponding  to  different  events  are 

p^  =  Pr{event  of  type  j  |  transition  (occurrence  of  Poisson  (A) 

process)  when  in  state  s)  =  r~*/A  j  =  1,2,...,E 

s 

p^  =  Pr{null  event,  i.e.,  remaining  in  state  s  |  transition  when 

in  state  s}  =  1  -  r  /A 

s 

t^  =  the  state  to  which  system  goes  when  a  null  event  occurs  -  s  . 


and 
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Transient  probabilities  of  the  continuous-time  Markov  process  are  com¬ 
puted  from  the  transient  probabilities  of  the  Markov  chain,  using  equa¬ 
tion  (8), 

T(e , t) 

tt( t)  =  £  4>(n )  P(N(t)  =  n)  . 

n=0 


The  vectors  4>(n)  are  computed  recursively: 


(i)  <j>(0)  =  tt(0) 

(ii)  cj)(n+l)  is  computed  from  cj>(n)  as  follows: 

(a)  (J)  (n+1)  =  <j>  (n)  •  p^,  s  e  S 

s  s  s 


then 


(10) 


(b)  for  j  =  1,2,..., G,  and  for  s  E  S,  add 

(j)  (n)  •  p-^  to  <j)  .(n+1)  . 

a  s 

s 

For  a  sequence  of  times  0  <  t^  <  <  . . .  <  over  which  the  process  is 

time-homogeneous,  compute  Tr(t  ),  m  =  1,...,£  simultaneously  by  initial- 

~  m 

izing  each  to  0  and  adding  d>(n)P(N(t  )  =  n)  to  the  mth  vector  as  the 

-  ~  m 

<J>(n)'s  are  computed  recursively.  At  a  discrete  point  of  time-nonhomogeneity, 
use  the  last  7r(t^)  as  the  initial  vector  and  repeat  the  process  for  time 
points  in  the  next  time  interval  for  which  time  homogeneity  exists. 

In  the  following  subsection,  we  illustrate  these  concepts  on 
three  examples,  the  first  dealing  with  a  classical  machine  repair  prob¬ 
lem,  the  second  an  M/E^/c  queue,  and  the  last  a  maintained  reliability 
system. 
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5.1 _ A  Machine  Repai  r  Example 


I 

example 

process 


Consider  the  classical  machine  repair  w:  th  spares  model  (see,  for 
,  G ross  and  Harris  [1974  }).  Since  this  is  a  finite  birth-death 
,  the  generator  is  given  by  equation  (1),  where 


!MA 

(N  -  s) A 
0 


y 


s 


(0  <  s  <  Y) 

(Y  <  s  <  Y+M  E  N) 
(s  >  N) 

(0  <  s  <  c) 

(c  <  s  <  N) 


and 


1/A  =  mean  time  to  failure  of  a  machine 
1/y  =  mean  time  to  repair  of  a  machine 
M  =  desired  number  of  machines  operating 
Y  =  number  of  spare  machines 
c  =  number  of  repairmen  (service  c'.ic  .eis)  . 

Let  us  consider  a  specific  problem  for  Wi ca  M  =  ?  =  3,  c  =  2, 

A  =  0.2,  and  y  =  0.4.  Thus, 


-1.0  1.0  0 

.4  -1.4  1.0 

0  .8  -1.8 

0  0.8 


0  0  0 

0  0  0 

1.0  0  0 

-1.8  1.0  0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 


Q  = 


0 


0  0 


.8  -1.6  .8  0 


0  0 


0  0  0  0 


.8  -1,4  .6  0 


0 


0  0  0 
0  0  0 
0  0  0 


0  0.8 

0  0  0 

0  0  0 


-1.2  .4  0 

.8  -1.0  .2 

0  .8  -.8 
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The  machine  repair  problem  has  a  special  structure  that  can  be 
used  to  develop  an  implementation  of  the  randomization  algorithm  that 
exploits  the  sparseness  of  Q,  which  is  common  to  many  applications  in 
queuing  theory:  there  is  a  very  small  class  of  distinct  events  that  are 
responsible  for  all  the  state  changes  of  the  system,  namely,  there  are 
only  two  types  of  events,  "failure  of  a  machine,"  and  "completion  of 
repair  of  a  machine."  If  the  state  of  the  system  is  described  as  the 
number  of  machines  in  repair,  say  s,  then  s  is  an  element  of  the  state 
space  S  (s  £  S),  where  S  =  {0,1, 2, . . . , 8}  in  our  example,  the  event 
"failure"  corresponds  to  transition  s  -*■  s  +  1  ,  the  event  "repair" 
corresponds  to  transition  s  -*  s  -  1  ,  and  the  respective  transition 
rates  (or  intensities)  are  found  on  the  superdiagonal  and  the  subdiag¬ 
onal  of  Q.  Thus  instead  of  describing  the  stochastic  behavior  of  this 
system  with  a  9  x  9  *  81  element  Q  matrix,  we  can  describe  it  in  terms 
of  the  state  space  S  =  {0,1,..., 8}  ,  two  events  E  =  {f,r}  — "failure" 
and  "repair" —  and  for  each  event  a  nine-component  vector  of  transition 
rates,  namely,  r^  and  rr,  and  a  nine-component  vector  indicating  the 
resulting  (target)  state  when  each  event  occurs  to  the  system  in  each 

•  f  IT 

possible  state,  namely,  t  and  t  .  Thus  the  vector  of  transition  rates 
with  which  a  "failure"  occurs  is 

rf  =  (1.0,  1.0,  1.0,  1.0,  .8,  .6,  .4,  .2,  0.0) 

and  the  vector  of  target  states  given  a  failure  is 

/  =  (1,  2,  3,  4,  5,  6,  7,  8,  8) 

(where  we  use  the  convention  that  the  target  is  the  same  as  the  original 

if  the  transition  intensity  is  0.0).  Similarly,  for  the  "repair"  event 
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the  transition  rates  are 


rr  =  (0.0,  .4,  .8,  .8,  .8,  .8,  .8,  .8,  .8) 
and  the  targets  are 

tr  =  (0,  0,  1,  2,  3,  4,  5,  6,  7)  . 

Thus  the  behavior  of  the  process  can  be  described  using  4  x  9  =  36 
numbers  instead  of  9  *  9  =  81  elements  of  Q.  This  description  of  the 
process  in  terms  of  the  state  space,  events,  transition  rate  vectors, 
and  target  vectors  is  our  SERT  procedure. 

We  implement  the  randomization  algorithm  using  this  idea  of  rate 
vectors  and  target  vectors  for  each  event,  cf.  Grassmann's  QUE  package 
(Grassmann  [1982]).  Consider  the  Markov  chain  with  transition  probabil¬ 
ity  matrix  P  on  the  underlying  Poisson  process:  it  has  three  events, 
"failure,"  "repair,"  and  "null  event"  (no  transition) .  We  can  see  that 
this  Markov  chain  can  be  described  in  terms  of  transition  probability 
vectors  for  each  event  and  vectors  of  target  states  for  each  event. 

Note  that  the  vectors  of  targets  are,  of  course,  the  same  for  the  orig¬ 
inal  Markov  process  and  the  Markov  chain  on  the  underlying  Poisson  pro¬ 
cess.  The  vectors  of  transition  probabilities  for  the  three  events  are: 


t 

p  = 

(.556, 

.556,  .556,  .556, 

.444, 

.333, 

.222, 

.111, 

0) 

r 

P  = 

(0, 

.222,  .444,  .444, 

.444, 

.444, 

.444, 

.444, 

.444) 

(.444, 

.222,  0,  0, 

.112, 

.222, 

.333, 

.444, 

.556)  . 

Note  that  the  vectors  p^,  pr,  p^,  t^,  tr,  and  the  value  A  =  1.8  completely 

f  f  r  r 

determine  the  process.  Note  that  p  =  r  /A  and  p  =  r  /A  ,  where 
r^  is  the  superdiagonal  and  rr  is  the  subdiagonal  of  the  Q  matrix. 
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The  fundamental  relationship  for  the  randomization  algorithm  to 
be  used  is  equation  (8),  namely, 

T 

7T(t)  =  £  <J>(n)P(N(t)  =  n)  . 

i=0  ~ 

The  vectors  4>(n) ,  n  =  1,2,...,  are  computed  recursively;  <{>  =  ^(0), 

which  is  assumed  given.  To  compute  <Ji(n+l)  from  4>(n)  ,  recall  that 
<J>(n+l)  =  <J)(n)P.  We  want  to  avoid  direct  matrix  multiplication  by  P 
because  it  is  so  sparse.  Instead  we  use  the  algorithm  given  as  (10) 
based  on  the  vectors  of  transition  probabilities  p^,  pr,  p^,  and  the 
vectors  of  target  states  t  and  t  . 

We  restate  the  algorithm  for  this  example  as  follows.  For  this 
system  there  are  three  steps: 

1.  For  each  s  e  S,  set  <J>  (n+1)  =  <J>  (n)  •  p^ 

s  s  s 

2.  For  each  s  e  S,  add  <J)  (n)  •  p^  to  <J>  ,(n+l) 

S  S 

s 

■J“ 

3.  For  each  s  e  S,  add  <j>g(n)  *  Pg  to  <{>  r(n+l) 

The  notation  takes  a  bit  of  getting  used  to.  Subscripts  always  refer  to 
an  element  of  a  vector.  For  example, 

4>g (n)  =  sth  element  of  <J>(n) 

=  Pr{system  state  is  s  (s  machines  in  repair)  after  n 
occurrences  of  the  underlying  Poisson  (A)  process) 

t^  =  sth  element  of  the  failure  target  vector,  t^;  that  is, 

the  state  to  which  the  system  transits  given  it  is  in 
state  s  and  a  failure  occurs 
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<p  c(n+l)  =  t^fh  element  of  the  <J>(n+l)  vector 

t  s  ; 

=  Prlsystem  state  is  t  after  n  +  1  occurrences  of  the 
underlying  Poisson  (A)  process}  . 

A  little  thought  reveals  that  this  algorithm  is  exactly  the  same 
as  <J>(n)P,  but  avoids  all  the  multiplication  by  zeroes  off  the  three 
main  diagonals,  and  also  avoids  storage  of  an  N  x  N  P  matrix. 

Consider  the  P  matrix  for  our  example.  We  can  calculate  <J)(n+l)  = 
(n) P  by  straight  matrix  multiplication.  For  example, 

<t>4  (n+1)  =  0*4»0(n)  +  O'^Cn)  +  0*(J>2(n)  +  .556<J>3(n)  +  ,112<f>^(n) 

+  .444<t>,-(n)  +  0* <{>,  (n)  +  0* cf>^(n)  +  0*<J>_(n) 

=  .556<J)^(n)  +  .  1124>^ (n)  +  .444<J),_(n)  . 

Notice  all  the  zero  multiplications.  The  algorithm  method  builds  up  the 
<j>(n+l)  vector  as  follows: 


Algorithm 

Step 

<j>(n+l)  Element 

(J)Q(n+l) 

^(n+1) 

♦2(n+l) 

4>3(n+l) 

(J)^  ( • 

1 

.444<J>Q(n) 

.222^(11) 

0*<f>2(n) 

0*<{)3(n) 

.  1124>^  (n) 

+ 

+ 

+ 

+ 

+ 

2 

.5564>0(n) 

.556^(0) 

.556<j)2(n) 

.5564>3(n)  ... 

+ 

+ 

+ 

+ 

3 

0+.222<{)1(n) 

•  444<J>2(n) 

.4444>3(n) 

.444<)>4  (n) 

.4444>3(n) 

The  column  sums  directly  give  <t>g(n+l);  for  example, 

<(>4(n+l)  -  .112<j>4(n)  +  .556<(>3(n)  +  .444<j>5(n) 
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and  hence  the  six  zero  multiplications  are  avoided. 

For  the  machine  repair  problem,  this  algorithmic  approach  has 
been  written  in  a  FORTRAN  program  called  REPTRAN1  (see  Gross,  et  at. 
[1982]).  It  can  handle  discrete  changes  in  failure  and  repair  rates. 
Over  an  interval  for  which  Q  is  constant  it  simultaneously  computes 
Tr(t)  for  a  lattice  of  time  points,  by  initializing  the  vectors  to  0  and 
then  adding  <Ji(n)P(N(t)  =  n)  as  the  4>(n)'s  are  computed  recursively, 
until  the  truncation  point  is  reached.  The  transient  probability  vec¬ 
tors  for  the  next  interval  of  constant  Q  is  computed  using  the  prob¬ 
ability  vector  of  the  last  point  in  the  preceding  interval  as  the  initial 
vector  for  this  interval.  (It  computes  A  from  the  formula 

SMX  +  cy  if  c  <  Y 

MX  +  Yy  if  c  >  Y,  X  >  U 

(M  +  Y  -  c)A  +  cy  if  c  >  Y,  A  <  y 

instead  of  searching  for  the  maximum  value  of  qg.)  It  also  plots  avail¬ 
ability  (the  probability  of  at  least  M  machines  up  =  1  -  probability  of 
at  most  Y  down)  versus  time.  Figure  1  shows  the  plot  for  this  example, 
where  X  increases  by  50%  at  time  6  and  y  increases  by  50%  at  time  10. 

5.2  A  Queuing  Example;  M/ E c  System 

In  this  section  we  give  a  SERT  model  of  an  M/E^/c  queuing  system. 
The  transient  state  probabilities  can  then  be  solved  using  the  randomi¬ 
zation  algorithm  given  in  the  previous  section. 

First  consider  the  state  description  of  the  system.  One  possi¬ 
bility  is 

s  —  (s^ ,  • * • »  q) 
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Figure  1. — Availability  vs.  time  for  machine  repair  example  with  a  50% 
increase  in  \  at  t  =  6  and  a  50%  increase  in  y  at  t  *  10. 
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where  equals  the  Erlang  stage  of  the  1th  service  channel  (s^  =  0 
signifying  that  the  1th  channel  is  empty)  and  q  equals  the  number  of 
customers  waiting  in  the  queue.  We  can  partition  the  state  space  into 
two  sets  corresponding  to  whether  the  queue  is  empty,  S^,  or  whether 
there  is  a  positive  number  of  customers  in  the  queue,  S+.  For  states  in 
Sq,  q  must  be  0  but  each  s^  may  be  any  integer  between  0  and  k;  thus 
consists  of  (k  +  1)  states,  which  we  assume  to  be  arranged  in  lexico¬ 
graphical  order.  Since  the  queuing  system  has  infinite  capacity,  S+  is 
clearly  infinite  and  thus  must  be  truncated  in  order  to  use  the  randomi¬ 
zation  algorithm,  as  discussed  in  Section  4.6.  Let  S  correspond  to 

+,K. 

the  subset  of  S  for  which  q  _<  K  -  c  ;  thus  S  U  S  are  the  states 

*  0  TjK 

of  an  M/E, /c/K  system.  The  elements  in  S  must  have  1  <  q  <.  K  -  c  and 

K.  TjK  — 

1  <  s  <  k  ,  and  consequently  there  are  (K  -  c)k  states  in  S  ,  which 

1  i* y 1\ 

we  also  assume  are  arranged  in  lexicographical  order.  Finally,  we  shall 
add  one  more  state,  an  absorbing  state  "a",  which  will  be  useful  in  ob¬ 
taining  a  bound  on  the  error  due  to  truncation.  Let  S  =  {a}.  Whenever 

a 

the  system  exceeds  K  customers,  it  is  absorbed  into  S  .  It  will  be  nec- 

a 

essary  to  choose  K  by  observing  the  amount  of  probability  that  is  absorbed 

by  S  and  adjusting  K  to  make  this  equal  to  or  less  than  some  small 
3 

-2  -3 

value  such  as  10  or  10  .  The  state  space  will  be 

S-SK  =  S0VS+,KUSa- 

Now  let  us  consider  the  set  of  possible  events  E.  There  are  two 
basic  types  of  events:  arrivals  and  completions  of  a  stage  of  service. 
There  is  the  possibility  of  describing  several  different  systems  based 
upon  the  mechanism  by  which  a  customer  selects  a  service  channel  when 
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more  than  one  is  idle.  We  shall  model  this  as  a  classical  overflow  prob¬ 
lem:  the  arriving  customer  attempts  to  enter  channel  1;  if  channel  1  is 

busy  he  tries  channel  2,  and  so  forth.  If  all  channels  are  busy  he 
overflows  into  the  queue.  (In  our  truncated  state  space,  if  K  customers 
are  in  the  system,  the  arrival  forces  the  system  into  the  absorbing 

state  S  .)  Thus  the  state  change  due  to  an  arrival  is  deterministic 
a 

and  it  suffices  to  have  a  single  arrival  event  A.  It  is  necessary  to 
have  a  separate  event  of  completion  of  a  stage  of  service  for  each  ser¬ 
vice  channel,  because  the  target  state  depends  on  the  channel  where  the 
event  occurred.  Thus  we  have  stage  completion  events  »...,Cc  and 

the  event  set  has  c  +  1  elements 

E  —  (a,  ,  ^2*  *  *  *  *  }  . 

Now  let  us  consider  the  collection  R  of  rate  vectors.  We  have 
c  +  1  rate  vectors,  one  corresponding  to  each  event  in  E.  In  the  model¬ 
ing  activity  we  think  in  terms  of  a  rate  function,  i.e.,  the  transition 
rate  is  a  function  of  the  state  of  the  system.  (In  Sections  5.2  and  5.3 

I 

we  thus  use  functional  notation  instead  of  vector  notation  to  describe 

the  rates  and  targets  in  the  SERT  description.)  First  consider  the  rate 

function  for  the  arrival  event  A:  arrivals  keep  coming  regardless  of 

the  state  of  the  system;  thus 

r.(s)  =  A  s  e  S  . 

A 

In  our  standard  numerical  implementation  of  the  randomization  procedure 
we  would  use  the  above  function  to  generate  the  rate  vector  with  compo¬ 
nents  rg  *  ^  *  s  e  5,  i.e., 

rA  =  (A,A,...,A) 
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which  would  be  used  in  the  calculation.  (There  is  always  some  flexi¬ 
bility  in  specifying  the  SERT  description.  In  this  case  we  could  have 
set  r  (a)  =  0  instead  of  X,  but  it  will  make  the  programming  easier 

A 

to  have  the  arrival  rate  equal  to  X,  independent  of  the  state 
of  the  system.  A  similar  consideration  holds  for  the  rate  of  service- 
stage  completion.)  Now  consider  the  event  C^,  completion  of  a  stage  of 
service  at  channel  i.  Since  the  server  in  channel  i  completes  stages 
according  to  a  Poisson  process  with  rate  y^, 

rc  (•)  -  ut  s  e  S  . 

(Again  it  is  convenient  for  programming  purposes  to  have  a  constant  rate. 
In  this  case  if  the  server  completes  a  stage  of  service  while  the  chan¬ 
nel  is  empty,  it  remains  empty,  constituting  a  dummy  event.) 

Finally,  let  us  consider  the  collection  of  all  target  vectors  or 
functions,  T.  Each  of  these  gives  the  target  state  as  a  function  of 
current  state  for  each  event  in  E.  First  consider  the  arrival  event. 
Since  we  are  modeling  the  system  as  a  classic  overflow  system 
t^(s^,S2» . .  .  ,sc»  q)  —  (Sj^, . . .  ,s^  ^ )  1,  s  ^  |  ^  * •  •  • » 

if  s^  >  0,  j  =  l,...,i-l,  and  s^  =  0 

=  (s. ,...,s  ,  q  +  1) 

I  c 

if  s^  >  0,  j  =  l,...,c,  and  q  <  K  -  c 

=  a  if  q  =  K  -  c 

tA(a)  =  a  . 

Now  consider  the  service-stage  completion  events.  Whenever  this  event 
occurs  the  channels  move  to  the  next  stage  unless  it  is  idle 
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t  (s1 ,  >  .  ■  ,s  ,  q  )  -  ( s  1  ,  .  .  .  ,  s  .  "fl 1  •  •  m s  i  q  ) 

1  c  11  c 

if  0  <  s .  <  k 

l 

=  (s^  * • • •»  ^  * • *  * » (  q  —  1) 
if  s .  =  k  and  q  >  0 

l 

=  (s^, . . . ,0, . . . ,s^,  q) 

if  s.  =  k  and  q  =  0,  or  if  s,  =  0 
l  i 

tc/a)  -  a 
i 

This  completes  the  SERT  description  of  the  model.  It  is  a  straight¬ 
forward  programming  task  to  construct  the  rate  and  target  vectors  corres¬ 
ponding  to  the  lexicographically  ordered  state  space  and  to  compute  tran¬ 
sient  state  probabilities  using  the  randomization  algorithm.  There  will 
be  two  sources  of  probability  loss:  the  e  probability  lost  due  to  trunca¬ 
tion  of  the  infinite  series  (4),  and  the  probability  absorbed  in  a  , 
it  (t)  .  The  total  error  on  all  probabilities  computed  will  be  bounded  by 

Si 

-(£  +  it  (t))  and  0.  If  this  bound  is  not  satisfactory,  £  can  be  decreased 

Si 

and  K  increased. 


5.3  A  Reliability  Example:  Parallel 
System  with  Component  Repair 

Consider  a  system  of  n  components  in  a  parallel  configuration. 
Initially  all  components  are  operating.  The  components  have  independent 


exponentially  distributed  lifetimes;  the  failure  rate  of  the  ith  component 


is  Upon  failure  the  component  is  repaired.  The  repair  time  has  an 

independent  generalized  Erlang  distribution  with  k.  stages  and  rate  y.  . 

i  l » j 

of  completion  of  the  ith  stage  of  repair,  j  =  l,2,...,k^.  After  repair 
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is  completed  the  component  functions  for  another  independent  exponential 
(A^)  period  of  time  before  failing  again.  System  failure  occurs  when  all 
n  components  are  simultaneously  under  repair.  Let  L  equal  the  time  until 
system  failure.  We  would  like  to  compute  Pr{L  <  t)  .  Analytical  and 
approximate  aspects  of  this  and  related  problems  are  given  by  Brown  [1975] 
and  Keilson  [1975]. 

First  let  us  consider  the  ith  component.  This  component  can  be  in 
any  of  k^  +  1  states  (s^  =  0,1,2, ... ,k^) ,  where  state  0  corresponds  to 
a  functioning  unit  and  j  corresponds  to  the  jt?z  stage  of  repair, 
j  =  l,2,...,k^  .  The  Markov  process  describing  the  ith  component  has 
generator 


This  process  has  one  basic  type  of  transition:  "change."  Thus  we  will 
be  able  to  describe  the  evolution  of  the  n-component  system  with  n 
events:  E  =  {C^,  C2,  .  ..,  Cn)  ,  where  is  the  event  of  "change"  of 

the  ith  component. 

Now  let  us  consider  the  state  of  the  entire  system.  Let 

s  =  (s^^, . .  •  ,sn) 

where  s^  equals  the  state  of  the  ith  component,  s^  =  0,l,2,...,k^  , 
for  i  =  1, . . . ,n.  Let 
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S*  =  {(s,,...,s  )  :  0  <  s.  <  k.,  i  =  l,...,n) 
1  n  =11 

5  —  {(s  , . . . ,s  )  ;  1  ^  s .  5  k  • »  ^  = 

r  1  n  1  1 


S  =  S*  _  S 
0  *  ‘ 


The  subset  S  is  the  set  of  all  states  corresponding  to  system  failure 

r 

and  the  subset  is  the  set  of  all  states  corresponding  to  an  operating 
system.  Assume  that  is  ordered  lexicographically.  Let  {f}  be  a 
single  absorbing  state  corresponding  to  system  failure.  The  state  space 
for  this  process  can  be  described  as 


5  =  SQ  V  {f }  . 

Now  let  us  consider  the  collection  R  of  rate  vectors,  one  for  each 
event  type  in  E.  Consider  the  event  The  rate  vector  (or  function) 

for  this  event  is 


rC1(sl,S2’,**’Sn)  “  Ai 


if  s±  =  0 


rr  (s, *sn, . • • ,S  )  =  n.  if  s  >  0 

.  x  z  n  i  *  s .  i 

l  *  l 


1 


Finally,  the  target  vectors  are 


tr  (s .,..., s  )  (s1 , . . . ,s .+1, . . . ,s  ) 

x  n  x  i  n 

if  0  <  <  k  or 


if  s,  <  k  and  s^  =  0  for  some  i  ?  j 


(s^, . . . ,0, . . . ,sn) 


if  s.  =  k, 
i  i 


if  s  =  0  and  s.  >0  for  all  i  ^  i 
i  J  J 
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tc  (f)  =  f 
1 

This  completes  the  SERT  description  of  the  process.  It  can  be  programmed 
and  the  randomization  technique  used  to  compute  ^(t),  which  will  equal 
Pr{L  <  t}  . 

6.  COMPUTATION  OF  OTHER  TRANSIENT  QUANTITIES 
Using  the  randomization  construction,  useful  computational  formu¬ 
lae  can  be  developed  for  other  transient  quantities:  (i)  expected 
time  averages,  (ii)  first  passage  time  distributions,  and  (iii)  expected 
number  of  events  of  a  certain  type  occurring  during  a  time  interval. 

b.l  Expected  Time  Averages 

Suppose  {X(t),  t  >  «>}  is  a  continuous  time  Markov  process  with 
generator  Q,  where  X(t),  the  state  that  the  system  is  in  at  time  t,  can 
take  values  s  =  0,1,2,...,N  .  Further  suppose  that  A  =  max|qg|  <  00  , 
i.e.,  the  process  is  uniformizable .  Then  from  the  preceding  discussion, 
the  process  can  be  represented  as  a  discrete  time  Markov  chain  on  an 
underlying  Poisson  process,  with  rate  A  and  transition  matrix  P  = 

Q/A  +  I  .  Let  (N(t),  t  >_  0}  be  the  counting  process  of  the  underlying 
Poisson  process,  i.e., 

N(t)  =  it  occurrences  of  Poisson  process  in  [0,t] 

and  let  {T  ,  n  =  1,2,...}  be  the  times  of  occurrence  of  this  Poisson 
n 

process.  (Then  T  =  min(t  :  N(t)  >  n)  ,  etc.)  The  discrete  time  pro- 
n 

cess  {Y  ,  n  =  0,1,2,...}  ,  then,  is  the  state  of  the  system  just  after 
n 

the  nth  occurrence  of  the  Poisson  (A)  process.  Thus  (Y  ,  n  =  0,1,...} 

n 
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is  the  discrete  time  Markov  chain  with  transition  matrix  P.  With  this 
terminology  we  can  now  proceed  to  the  quantities  of  interest. 


00  n 


THEOREM 


i.  E[jC5  X(T)dT]  =  l  l  EY.  /  (n  +  1)  (e  At(At)n)/n!  . 

n=0  lj=0  3 


Proof :  E[^)t  X(x)dx] 


=  l  E[jft  X(x) dx  |  N(t)  -  n]  P (N(t)  =  n) 
n=0  U 

OO 

=  l  EIYq^  +  Y2(T2-T  )  p  ...  +  *  (t-T  )  I  N(t)  =  n]  P(N(t)  =  n 
n=0 

=  I  e|  I  Y.(T  -  T  )  |  N(t)  =  n  |  P(N(t)  =  n) 
n=0  Lj=0  3  3  3  J 

°°  n 

=  11  E[y j <Tj+1  -  Tj)  |  N(t)  =  n]  P (N(t)  =  n) 


n=0  j=0 


oo  n 

=  l  l  EY,  E[T  -  T  I  N(t)  =  n]  P (N(t)  =  n) 
n=0  j=0  J  J  1  J 

a,  n 

=  I  I  EY.  1/ (n+1)  P (N(t)  =  n) 
n=0  j=0  J 

00  n  . 

-  l  l  EY.  /  (n  +  1}  (e"At(At)n)/n! 
n=0  ij=0 


Q.E.D 


The  associated  computational  formula  in  terms  of  the  transient 
probability,  <p,  for  the  Markov  chain  will  be 


oo  n 


E[^5  X(x)dx]  =  l  l  EY.  (e"At(At)n)/(n  +  1) ! 
n=0  j=0  3 


=  l  ll  s<f>  As)  (e‘At(At)n)/(n  +  1)!  . 
n=0  I j=0  seS  3 
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The  above  results  can  be  generalized  in  the  obvious  straightfor¬ 
ward  way  to  computational  formulae  for 

El  £  f(X(x)>dx] 

where  f  is  any  real  valued  function  on  the  state  space  of  the  process. 

As  an  example,  consider  a  machine  repair  problem.  Let  X(x)  equal 
the  number  of  good  machines  available  at  time  T.  Let  f (x)  =  min(x,M)  , 
where  M  is  the  desired  number  of  operating  machines.  Then 

El/J  min(X(x),M)dx] 

is  the  expected  amount  of  machine  operating  time  achieved  during  [0,t]. 

As  a  similar  example,  consider  a  fleet  of  aircraft  of  which  M  will  be 
kept  in  flight  operations  if  M  are  available.  Then  the  integral  repre¬ 
sents  the  amount  of  flight  operation  time  achieved  by  the  fleet  during 
[ 0 , t 1 .  If  Ts  equals  the  sum  of  deterministic  flight  time  per  sortie  and 
the  deterministic  inter-sortie  ground  time,  then 

Efjf  min(X(x),M)dx]  /  Tg 

equals  the  expected  number  of  sorties  flown  in  [ 0 , t ] .  (This  assumes  that 
the  denominator  is  much  less  than  the  numerator,  otherwise  it  is  only  an 
approximation  unless  t  is  an  integer  multiple  of  sortie  time,  in  which 
case  it  is  always  exact.) 

Other  performance  measures  that  can  be  computed  in  a  similar  way 
include  the  expected  time  the  process  spends  in  a  particular  subset  of 
the  state  space,  say  A,  during  [ 0 , t ] : 

Eljf  lA(X(x))dx] 
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where  l.(*)  is  the  indicator  function  of  the  set  A.  If  A  indicates  the 
A 

set  of  states  for  which  an  availability  requirement  is  satisfied,  then 
the  integral  represents  the  amount  of  time  during  [0,t]  that  the  require¬ 
ment  is  met . 

6.2  First  Passage  Time  Distributions 

Suppose  (X(t),  t  >  0)  is  a  Markov  process  on  the  finite  state 
space  S.  Let  A  be  a  subset  of  S  and  TA  the  first  passage  time, 

T  =  min{t  :  X(t)  £  A}. 

A 

The  randomization  technique  can  be  used  to  compute  passage  times  as  fol¬ 
lows:  define  an  associated  process  {X  (t) ,  t  ^  0}  on  the  state  space 

Pi 

S  -  A  \J  {a}  with  generator  Q  such  that  q  .  .  =  q.  .  for  i,j  t  A, 

A  A;i, j  1 » J 

q.  .  =  \  q.  .  for  i  t  A,  and  q  .  =  0  for  all  j.  Thus  the  set  A 

’  ,a  jf.A  A;a»J 

is  collapsed  down  to  one  absorbing  state  a  .  The  transition  intensities 

are  identical  except  for  0  rate  out  of  A.  Thus 

P{T  <  t}  =  P(X. (t)  =  a) 

A  —  A 

and  the  first  passage  probability  (left-hand  side  of  the  equation)  is 
equivalent  to  computing  a  transient  state  probability  for  a  Markov  pro¬ 
cess  (generator  Q  ) ,  which  can  be  computed  using  the  randomization  algo- 
rithm.  Melamed  and  Yadin  [1980,1981]  use  this  approach  to  compute  so¬ 
journ  times  for  queuing  networks. 

6.3  Expected  Number  of  Events  in  an  Interval 


example,  in  a  reliability  problem  we  might  be  interested  in  the  expected 
number  of  failures  occurring  in  [0,t]. 


Let  e^  £  E  and  suppose  we  are  computing  $(n)'s  recursively, 

where  <}>  (n)  =  P(X(T  )  =  s)  =  P(Y  =  s)  .  Then  <p  (n)  •  p,(s)  = 
s  n  n  s  j 

P(Y  =  s  and  next  event  is  e.)  .  Thus  the  inner  product  £  <f)  (n)  ' 


Pj(s)  equals  PCe^  occurs  at  Tn+^)  ,  which  also  equals  the  expected 

number  of  e.'s  in  (T  ,  T  , .  ]  .  Thus 
j  n  n+1 


E[ number  of  ej's  in  [ 0 , t ] ] 

00 

=  £  E[number  of  e.’s  in  [0,t]  |  N(t)  =  n]  P[N(t)  =  n] 

n=0  J 

00 

=  J  E[number  of  e.'s  in  [0,T  ]  |  N(t)  =  n]  P[N(t)  =  n] 
n=0  J  n 

oo  n-l 

=  J  [  E [ number  of  e.'s  in  (T,  ,T,  .]  |  N(t)  =  n]  P[N(t)  =  n] 
n=0  k=0  J  k  K  1 

00  n-l 

=  I  l  l  4>  (k)  •  P  (s)  P[N(t)  =  n]  . 
n=0  k=0  scS  ^ 


7.  A  GENERAL  IMPLEMENTATION  FOR  PROCESSES 
WITH  SPARSE  GENERATORS 

The  SB RT  approach  works  well  if  each  event  can  occur  while  the 
system  is  in  an  arbitrary  state  or  if  there  are  only  a  few  states  where 
a  given  event  cannot  occur.  For  example,  in  the  machine  repair  system, 
a  repair  cannot  occur  when  all  machines  are  functional  and  a  failure 
cannot  occur  when  all  machines  are  in  repair;  otherwise  each  event  can 
occur  in  all  remaining  states,  and  SERT  is  quite  efficient.  However, 
there  are  sparse  systems  that  are  not  amenable  to  this  approach,  e.g.. 


T 
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k  =  1 , 2 , . . . , n j  ,  j  =  1,2, ... ,N  .  All  the  information  in  a  generator  Q 
will  be  captured  in  the  vectors  r*  and  t*.  For  the  generator  in  equa¬ 
tion  (11)  this  representation  is  given  by 
N  *  5, 

A  -  20, 

N 

=  1,  n2  =  2,  n3  =  4,  n^  =  2,  n5  =  1,  N  +  J  ^  *  15, 

T1  =  {2},  T2  =  {3,4},  T3  =  (1,2, 4, 5},  T 4  =  {2,3},  =  {4}, 

t*  =  (1,  2,  2,  3,  4,  4,  1,  2,  4,  5,  2,  2,  3,  1,  4), 
r*  =  (-1,  1,  -6,  3,  3,  -20,  4,  4,  6,  6,  -12,  5,  7,  -2,  2), 
p*  =  (.95,  .05,  .7,  .15,  .15,  0,  .2,  .2,  .3,  .3,  .4,  .25,  .35,  .9,  .1)  . 
This  provides  an  economical  way  of  storing  a  sparse  generator  Q  of  mod¬ 
erate  or  large  size. 

The  transient  distributions  <£(n)  of  the  Markov  chain  on  the  under- 
lying  Poisson  must  be  calculated.  As  before,  this  is  done  recursively 
using 

4>(n  +  1)  =  <j>(n)P  .  (12) 

The  representation  of  Q  with  t*  and  r*  yields  a  representation  of  P  in 
terms  of  t*  and  p*  that  makes  recursive  calculation  of  4>  very  easy. 
Letting  (n  +  1)  be  PHINEW(J),  (n)  be  PHIOLD(J),  t*  be  TSTAR(l) ,  and 

p*  be  PSTAR(I),  a  FORTRAN  routine  for  equation  (12)  is: 
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I 

00  1  J  =  1,N 

|  1  PHINEW(J)  =  0.0 

I  =  0 

00  3  J  =  1,N 

|  1  =  1  +  1 

PHIJ  =  PHIOLD(J) 

PHINEW(J)  =  PHINEW(J)  +  PHIJ  *  PSTAR(I) 

|  MJ  =  TSTAR(I) 

IF  MJ.EQ.O  GO  TO  3 
DO  2  K  =  1,MJ 
1  =  1  +  1 
TJK  =  TSTAR(I) 

2  PHINEW(TJK)  =  PHINEW(TJK)  +  PHIJ  *  PSTAR(I) 

3  CONTINUE 

This  representation  and  algorithm  totally  exploit  any  sparseness  in  Q. 

For  large  systems  the  easier  approach  may  often  be  a  SERT  descrip¬ 
tion  of  the  process.  If  a  large  degree  of  sparseness  remains  in  the  SERT 
description,  i.e.,  many  rates  are  0,  the  set  of  target  vectors  T  can  be 
transformed  into  one  vector  t*  and  the  set  of  rate  vectors  R  can  be 
transformed  into  one  vector  r*  and  the  general  approach  used. 

In  some  instances  a  hybrid  approach  might  be  best.  If  some  events 
in  the  SERT  description  can  occur  at  all  states,  these  events  can  be 
modeled  with  individual  rate  and  target  vectors.  If  other  events  in  the 
SERT  description  are  impossible  at  many  states,  then  these  events  should 
be  pooled  into  one  general  "miscellaneous"  event  described  by  t*  and  r*. 

-  42  - 


In  closing,  we  note  that  we  are  concerned  with  systems  in  which 


a  large  number  of  occurrences  of  the  Poisson  process  will  happen  during 
the  transient  time  interval.  Thus  our  primary  concern  is  the  rapid  cal¬ 
culation  of  equation  (12) .  We  are  willing  to  spend  a  little  extra  time 
in  setting  up  a  one-dimensional  state  space  and  rate  and  target  vectors 
in  order  to  achieve  maximum  speed  in  this  calculation. 


8.  CONCLUSIONS  AND  SUMMARY 

We  have  presented  a  review  of  and  some  extensions  to  the  randomi¬ 
zation  technique  for  computing  transient  quantities  for  Markov  processes. 
We  have  introduced  the  general  modeling  approach  of  SERT  which  seems  to 
have  rather  broad  applicability  and  is  a  general  and  efficient  means  of 
implementing  the  randomization  technique.  Gross  and  Miller  (1982)  have 
used  the  SERT/Randomization  approach  on  multi-echelon  repairable  item 
provisioning  models.  Miller  (1982)  has  used  SERT  and  the  sparseness 
approach  of  Section  7  in  reliability  calculations  for  fault-tolerant 
reconf igurable  computer  systems. 

The  SERT  approach  has  the  advantage  that  it  is  general  and  applies 
to  many  systems.  However,  implementations  of  the  randomization  technique 
tailored  to  a  specific  system  may  be  more  efficient.  In  the  latter  two 
examples  of  Section  5  it  is  clear  that  various  improvements  could  be  made 
in  the  implementation  by  treating  special  cases.  For  example,  in  the 
M/E^/c  queue  it  is  unnecessary  to  store  entire  rate  vectors  because  they 
are  all  constant.  Grassmann  (1982)  has  written  an  efficient  package  to 
compute  state  probabilities  for  queuing  networks.  Melamed  and  Yadin 


!A 


(1981)  have  developed  an  efficient  algorithm  for  sojourn  time  distribu¬ 
tions  in  queuing  networks.  In  spite  of  this  tradeoff  between  general 
and  specific  algorithms,  it  appears  that  SERT  is  a  good  starting  point 
for  the  algorithmic  analysis  of  transient  Markov  processes. 
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